###setup
using Distributed
addprocs(4)

@everywhere using SharedArrays
@everywhere using Interpolations, Parameters, Distributions, DataFrames, CSV, LinearAlgebra, Statistics, Random, StatsBase, NLopt, Optim, Statistics

@everywhere dir = "C:\\Users\\Garrett\\Documents\\Grad_School\\Papers\\SIUM_JPE\\Replication\\Model"
@everywhere cd(dir)

@everywhere include("sium_model.jl")
@everywhere include("sium_utilities.jl")
@everywhere include("sium_estimation.jl")
@everywhere include("sium_simulations.jl")
@everywhere include("sium_readin.jl")
@everywhere include("sium_simulations_decomp.jl")

include("sium_make_weight_matrix.jl")
include("sium_make_marriage_plots.jl") 

####read in model utilities and targeted moments, and initialize model primitives
package = Readin_utilities(dir)
package_moments = Readin_moments(dir)

############Estimation Procedure
function wrapper(guess::Array{Float64}, grad::Vector)
    err = Objective_function(guess, package, package_moments; nsims = 20000, ret = 0, statcheck = 1)
    err
end

opt = Opt(:LN_SBPLX, 29)
opt.xtol_rel = 1e-4
opt.ftol_rel = 1e-4
opt.min_objective = wrapper

guess_init = [0.58673, 0.57223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]

err = Objective_function(guess_init , package, package_moments; nsims = 20000, ret = 0)
(minf, minx, ret) = NLopt.optimize(opt, guess_init)
println("guess_init = ", round.(minx, digits = 5))
println(minf)

err = Objective_function(minx, package, package_moments; nsims = 20000, ret = 0)
(minf2, minx2, ret) = NLopt.optimize(opt, guess_init)
println("guess_init = ", round.(minx2, digits = 5))
println(minf2)



########################Basic Simulation#######################
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]

output = Objective_function(guess_init, package, package_moments; nsims = 20000, ret = 1)

#@elapsed output = Objective_function(guess_init, package, package_moments; nsims = 20000, ret = 1, aopt = 4)
CSV.write("$dir\\Simulated_data\\simulated_data_base.csv", DataFrame(output[2], :auto), header=false) #write CSV file

########################Counterfactdual -- No Moving (w/o Behavioral Response)########################
#copy simulated data
data_clone = copy(output[2])
kid_pctiles = output[4]
data_clone[:,24] = data_clone[:,24]./data_clone[:,4]
data_clone[:,24] = data_clone[:,24].*data_clone[:,3] #recode

nrows = length(data_clone[:,1]) #useful for later
ncols = length(data_clone[1,:]) #useful for later

ranks = zeros(nrows, 1)
for i = 1:nrows
    ranks[i,1] = findmin(abs.(kid_pctiles .- data_clone[i,24]))[2]
end

data_clone[:,26] = ranks[:,1] #reset child percentile values
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_movecosts_nob.csv", DataFrame(data_clone, :auto), header=false) #write CSV file

########################Counterfactdual -- prohibitive moving costs########################
guess_cfact_move = copy(guess_init)
guess_cfact_move[10] = 750 #have to update this one
guess_cfact_move[19] = 750 #have to update this one too
prim_2000_cfact_move, res_2000_cfact_move = Initialize(package, package_moments)
prim_2010_cfact_move, res_2010_cfact_move = Initialize(package, package_moments; p2 = 1)
grids_cfact_move, est_cfact_move = Initialize_grids(prim_2000_cfact_move, guess_cfact_move) #initialize computing grids
V_iterate(prim_2000_cfact_move, est_cfact_move, res_2000_cfact_move, grids_cfact_move) #compute policy and value function.
V_iterate(prim_2010_cfact_move, est_cfact_move, res_2010_cfact_move, grids_cfact_move) #compute policy and value function.
data_simul_cfact_move = Simulate_cfact(prim_2000_cfact_move, prim_2010_cfact_move, est_cfact_move, res_2000_cfact_move, res_2010_cfact_move, grids_cfact_move, package, kid_pctiles; nsims = 20000) #simulate CHKS cohort, p25
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_movecosts.csv", DataFrame(data_simul_cfact_move, :auto), header=false) #write CSV file

guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]
prim_2000, res_2000 = Initialize(package, package_moments)
prim_2010, res_2010 = Initialize(package, package_moments; p2 = 1)
grids, est = Initialize_grids(prim_2000, guess_init) #initialize computing grids
V_iterate(prim_2000, est, res_2000, grids) #compute policy and value function.
V_iterate(prim_2010, est, res_2010, grids) #compute policy and value function.

data_simul_decomp_1 = Simulate_decomp_1(prim_2000, prim_2010, est, res_2000, res_2010, res_2000_cfact_move, res_2010_cfact_move, grids, package, kid_pctiles; nsims = 20000)
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_movecosts_decomp_1.csv", DataFrame(data_simul_decomp_1, :auto), header=false) #write CSV file

data_simul_decomp_2 = Simulate_decomp_2(prim_2000, prim_2010, est, res_2000, res_2010, res_2000_cfact_move, res_2010_cfact_move, grids, package, kid_pctiles; nsims = 20000)
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_movecosts_decomp_2.csv", DataFrame(data_simul_decomp_2, :auto), header=false) #write CSV file


###################################################################################################################
#######################Variance Decomp#######################
###################################################################################################################
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]

#demographic
package = Readin_utilities(dir)
package_moments = Readin_moments(dir)
package_decomp = Readin_decomp(dir)
package[1] = copy(package_decomp[1])
package[2] = copy(package_decomp[4])
package[4] = copy(package_decomp[7])
package[5] = copy(package_decomp[8])
package[6] = copy(package_decomp[9])
package[7] = copy(package_decomp[10])
package[8] = copy(package_decomp[11])
package[9] = copy(package_decomp[12])
package[10] = copy(package_decomp[13])
package[13] = copy(package_decomp[14])
output = Objective_function(guess_init, package, package_moments; nsims = 20000, ret = 1)
CSV.write("$dir\\Simulated_data\\simulated_data_decomp_d.csv", DataFrame(output[2], :auto), header=false) #


#economic
package = Readin_utilities(dir)
package_moments = Readin_moments(dir)
package_decomp = Readin_decomp(dir)
package[1] = copy(package_decomp[2])
package[2] = copy(package_decomp[5])
output = Objective_function(guess_init, package, package_moments; nsims = 20000, ret = 1)
CSV.write("$dir\\Simulated_data\\simulated_data_decomp_e.csv", DataFrame(output[2], :auto), header=false) #

#schooling
package = Readin_utilities(dir)
package_moments = Readin_moments(dir)
package_decomp = Readin_decomp(dir)
package[1] = copy(package_decomp[3])
package[2] = copy(package_decomp[6])
output = Objective_function(guess_init, package, package_moments; nsims = 20000, ret = 1)
CSV.write("$dir\\Simulated_data\\simulated_data_decomp_s.csv", DataFrame(output[2], :auto), header=false) #


###################################################################################################################
#######################Moving Cost Analysis########################
###################################################################################################################
package = Readin_utilities(dir)
package_moments = Readin_moments(dir)
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]
prim_2000, res_2000 = Initialize(package, package_moments)
prim_2010, res_2010 = Initialize(package, package_moments; p2 = 1)
grids, est= Initialize_grids(prim_2000, guess_init) #initialize computing grids
V_iterate(prim_2000, est, res_2000, grids) #compute policy and value function.
V_iterate(prim_2010, est, res_2010, grids) #compute policy and value function.
data_simul_movecosts = Simulate_movecosts(prim_2000, prim_2010, est, res_2000, res_2010, grids, package; nsims = 20000) #simulate CHKS cohort, p25
CSV.write("$dir\\Simulated_data\\simulated_data_movecosts.csv", DataFrame(data_simul_movecosts, :auto), header=false) #write CSV file

###################################################################################################################
#######################Counterfactual -- equalized school quality########################
###################################################################################################################
package = Readin_utilities(dir)
package_moments = Readin_moments(dir)
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]
output = Objective_function(guess_init, package, package_moments; nsims = 20000, ret = 1)
kid_pctiles = output[4]

####equalized school quality
package_cfact_g = Readin_utilities_cfact_1(dir)
package_moments = Readin_moments(dir)
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]
prim_2000_cfact_g, res_2000_cfact_g = Initialize(package_cfact_g, package_moments)
prim_2010_cfact_g, res_2010_cfact_g = Initialize(package_cfact_g, package_moments; p2 = 1)
grids_cfact_g, est_cfact_g = Initialize_grids(prim_2000_cfact_g, guess_init) #initialize computing grids
V_iterate(prim_2000_cfact_g, est_cfact_g, res_2000_cfact_g, grids_cfact_g) #compute policy and value function.
V_iterate(prim_2010_cfact_g, est_cfact_g, res_2010_cfact_g, grids_cfact_g) #compute policy and value function.
data_simul_cfact_g = Simulate_cfact(prim_2000_cfact_g, prim_2010_cfact_g, est_cfact_g, res_2000_cfact_g, res_2010_cfact_g, grids_cfact_g, package_cfact_g, kid_pctiles; nsims = 20000) #simulate CHKS cohort, p25
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_equal_g.csv", DataFrame(data_simul_cfact_g, :auto), header=false) #write CSV file


####equalized school quality, averaged
package_cfact_g_avg = Readin_utilities_cfact_3(dir)
package_moments = Readin_moments(dir)
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]
prim_2000_cfact_g, res_2000_cfact_g = Initialize(package_cfact_g_avg, package_moments)
prim_2010_cfact_g, res_2010_cfact_g = Initialize(package_cfact_g_avg, package_moments; p2 = 1)
grids_cfact_g, est_cfact_g = Initialize_grids(prim_2000_cfact_g, guess_init) #initialize computing grids
V_iterate(prim_2000_cfact_g, est_cfact_g, res_2000_cfact_g, grids_cfact_g) #compute policy and value function.
V_iterate(prim_2010_cfact_g, est_cfact_g, res_2010_cfact_g, grids_cfact_g) #compute policy and value function.
data_simul_cfact_g = Simulate_cfact(prim_2000_cfact_g, prim_2010_cfact_g, est_cfact_g, res_2000_cfact_g, res_2010_cfact_g, grids_cfact_g, package_cfact_g, kid_pctiles; nsims = 20000) #simulate CHKS cohort, p25
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_equal_g_avg.csv", DataFrame(data_simul_cfact_g, :auto), header=false) #write CSV file

#####equalized tuition rates
package_cfact_t = Readin_utilities_cfact_2(dir)
package_moments = Readin_moments(dir)
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]
prim_2000_cfact_t, res_2000_cfact_t = Initialize(package_cfact_t, package_moments)
prim_2010_cfact_t, res_2010_cfact_t = Initialize(package_cfact_t, package_moments; p2 = 1)
grids_cfact_t, est_cfact_t = Initialize_grids(prim_2000_cfact_t, guess_init) #initialize computing grids
V_iterate(prim_2000_cfact_t, est_cfact_t, res_2000_cfact_t, grids_cfact_t) #compute policy and value function.
V_iterate(prim_2010_cfact_t, est_cfact_t, res_2010_cfact_t, grids_cfact_t) #compute policy and value function.
data_simul_cfact_t = Simulate_cfact(prim_2000_cfact_t, prim_2010_cfact_t, est_cfact_t, res_2000_cfact_t, res_2010_cfact_t, grids_cfact_t, package_cfact_t, kid_pctiles; nsims = 20000) #simulate CHKS cohort, p25
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_equal_t.csv", DataFrame(data_simul_cfact_t, :auto), header=false) #write CSV file

#averaged
package_cfact_t_avg = Readin_utilities_cfact_4(dir)
package_moments = Readin_moments(dir)
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]
prim_2000_cfact_t, res_2000_cfact_t = Initialize(package_cfact_t_avg, package_moments)
prim_2010_cfact_t, res_2010_cfact_t = Initialize(package_cfact_t_avg, package_moments; p2 = 1)
grids_cfact_t, est_cfact_t = Initialize_grids(prim_2000_cfact_t, guess_init) #initialize computing grids
V_iterate(prim_2000_cfact_t, est_cfact_t, res_2000_cfact_t, grids_cfact_t) #compute policy and value function.
V_iterate(prim_2010_cfact_t, est_cfact_t, res_2010_cfact_t, grids_cfact_t) #compute policy and value function.
data_simul_cfact_t = Simulate_cfact(prim_2000_cfact_t, prim_2010_cfact_t, est_cfact_t, res_2000_cfact_t, res_2010_cfact_t, grids_cfact_t, package_cfact_t, kid_pctiles; nsims = 20000) #simulate CHKS cohort, p25
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_equal_t_avg.csv", DataFrame(data_simul_cfact_t, :auto), header=false) #write CSV file


#####loop for  shocks to g^ℓ
for i = 1:50
    println("working on ", i)
    @everywhere package = Readin_utilities(dir)
    @everywhere package_moments = Readin_moments(dir)
    package[1][i, 4] *= 1.1 #add 0.1 to 2000 skill price
    package[2][i, 4] *= 1.1 #add 0.1 to 2000 skill price
    package[1][i, 5] *= 0.9 #add 0.1 to 2000 skill price
    package[2][i, 5] *= 0.9 #add 0.1 to 2000 skill price
    @everywhere package_cfact_shock_g = copy(package)
    @everywhere guess_init = [0.59993, 0.507, -0.536, 0.375, 4.781, 0.885, 0.453, 0.31652, 0.28468, 11.771, 0.69359, 2.995, 1.07835, 2.00613, -1.714, 0.81437, 0.63844, 0.98475, 12.88188] ####FINAL FINAL FINAL
    output_simul = Objective_function(guess_init, package_cfact_shock_g, package_moments; nsims = 20000, ret = 1, short=1) #
    CSV.write("$dir\\Simulated_data\\simulated_data_cfact_shock_school_$i.csv", DataFrame(output_simul[2], :auto), header=false) #write CSV file
end

###################################################################################################################
########################Counterfactual -- skill price shocks########################
###################################################################################################################
#####loop for  shocks to g^ℓ
function simulate_shocks()
    for i = 1:50
        println("working on ", i)
        @eval @everywhere j = $i
        @everywhere package_temp = Readin_utilities(dir)
        @everywhere package_moments = Readin_moments(dir)
        @everywhere package_temp[1][j, 2] *= 1.1 #add 0.1 to 2000 skill price
        @everywhere package_temp[1][j, 3] *= 1.1 #add 0.1 to 2010 skill price
        @everywhere package_temp[2][j, 2] *= 1.1 #add 0.1 to 2000 skill price
        @everywhere package_temp[2][j, 3] *= 1.1 #add 0.1 to 2010 skill price
        @everywhere guess_init = [0.59993, 0.507, -0.536, 0.375, 4.781, 0.885, 0.453, 0.31652, 0.28468, 11.771, 0.69359, 2.995, 1.07835, 2.00613, -1.714, 0.81437, 0.63844, 0.98475, 12.88188] ####FINAL FINAL FINAL
        output_simul = Objective_function(guess_init, package_temp, package_moments; nsims = 20000, ret = 1, short = 1) #
        CSV.write("$dir/Simulated_data/simulated_data_cfact_shock_$i.csv", DataFrame(output_simul[2], :auto), header=false) #write CSV file
    end
end
simulate_shocks()


###################################################################################################################
########################Counterfactual -- Retention Policies########################
###################################################################################################################
function Simulate_retention(guess_init::Array{Float64,1}; nsims::Int64 = 20000)
    cfacts = [1 2 5]

    ####first go:

    #loop over states and counterfactual policies
    for l = 1:50, c in cfacts
        println("Working on state ", l, " and counterfactual policy ", c)

        #refresh
        @everywhere package = Readin_utilities(dir)
        @everywhere package_moments = Readin_moments(dir)
        output_simul = Objective_function(guess_init, package, package_moments; nsims = 20000, ret = 1, cfact = c, l_select = l, short = 1) #
        data_simul = output_simul[2]
        temp = length(data_simul[:,1])

        #add on counterfactual flags and push
        prim_2000, res_2000 = Initialize(package, package_moments)
        flag_cfact = ones(temp).*c
        flag_state = ones(temp).*prim_2000.fips[l]
        data_simul = hcat(data_simul, flag_cfact, flag_state)
        num = c * 1000 + l
        CSV.write("$dir\\Simulated_data\\simulated_data_cfact_retention_$num.csv", DataFrame(data_simul, :auto), header=false) #write CSV file
        finalize(res_2000)
        @everywhere GC.gc()
    end
end

guess_init = [0.59993, 0.507, -0.536, 0.375, 4.781, 0.885, 0.453, 0.31652, 0.28468, 11.771, 0.69359, 2.995, 1.07835, 2.00613, -1.714, 0.81437, 0.63844, 0.98475, 12.88188] ####FINAL FINAL FINAL
Simulate_retention(guess_init, nsims=20000)


###################################################################################################################
########################Identification Test########################
###################################################################################################################
function Identification_test(guess::Array{Float64,1})

    #re-initialize everything because why not
    package = Readin_utilities(dir)
    package_moments = Readin_moments(dir)
    guess_init = guess

    #baseline error
    error_base = Objective_function(guess_init, package, package_moments; nsims = 20000, ret = 0)
    #set up parameter and multiplier vectors
    modifiers = collect(-0.05:0.01:0.05)
    counter = 0
    data_id_test = Any[]

    for param in guess
        counter += 1
        #loop over multipliers
        for mod in modifiers

            #re-initialize parameter vector and modify appropriately.
            param_base = guess_init[counter]
            guess_test = copy(guess_init)
            guess_test[counter] = guess_test[counter]+mod
            error_test = 0.0

            #don't bother simulating if we're not changing anything
            if mod == 0.0
                error_test = error_base
            elseif mod!=1.0
                error_test = Objective_function(guess_test, package, package_moments; nsims = 20000, ret = 0)
            end

            #store in line: parameter we're looking at, value, error
            line = [param param_base+mod error_test]
            push!(data_id_test, line) #add to simulated data
        end
    end
    data_id_test = vcat(data_id_test...)
    data_id_test #return data
end

#obtain the data and store
guess_init = [0.59993, 0.507, -0.536, 0.375, 4.781, 0.885, 0.453, 0.31652, 0.28468, 11.771, 0.69359, 2.995, 1.07835, 2.00613, -1.714, 0.81437, 0.63844, 0.98475, 12.88188] ####FINAL FINAL FINAL
data_id_test = Identification_test(guess_init)
CSV.write("identification_test_data.csv", DataFrame(data_id_test, :auto), header=false) #write CSV file

###################################################################################################################
#functino for computing standard errors
###################################################################################################################
function standard_errors(guess::Vector{Float64}, perturb::Float64)
    num_sim = 20000 #10x normal simulation
    nrep = 1
    num_param = length(guess)
    num_moment = 234
    perturb_mat = zeros(num_moment, num_param)
    se_vec = zeros(num_param)

    #setup
    @everywhere package = Readin_utilities(dir)
    @everywhere package_moments = Readin_moments(dir)
    err_vec = Objective_function(guess, package, package_moments; nsims = num_sim, ret = 1)[5] #


    #loop over parameters
    for i = 1:num_param
        guess_perturb = copy(guess)
        guess_perturb[i] += perturb
        err_vec_perturb = Objective_function(guess_perturb, package, package_moments; nsims = num_sim, ret = 1)[5] #
        perturb_mat[:, i] = (err_vec.-err_vec_perturb)./perturb
    end

    gradient_mat = perturb_mat
    W = package[14]
    avar_mat_step1 = Matrix{Float64}(gradient_mat' * W * gradient_mat)
    avar_mat_step2 = pinv(avar_mat_step1)
    avar_mat_step3 = (1 + 1/nrep).*avar_mat_step2
    se_vec = sqrt.(diag(avar_mat_step3))
    se_vec
end

#get standard erors
guess_init = [0.59993, 0.507, -0.536, 0.375, 4.781, 0.885, 0.453, 0.31652, 0.28468, 11.771, 0.69359, 2.995, 1.07835, 2.00613, -1.714, 0.81437, 0.63844, 0.98475, 12.88188] ####FINAL FINAL FINAL
se_vec = standard_errors(guess_init, 1e-3)

for i = 1:length(guess_init)
    println(round(guess_init[i], digits=3), " || ", round(se_vec[i], digits=3))
end

############



###################################################################################################################
###################################################################################################################
###################################################################################################################
########################ROBUSTNESS SECTION########################
###################################################################################################################
###################################################################################################################
###################################################################################################################



###################################################################################################################
#######################Stationarity Assumption Check#######################
###################################################################################################################
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]
@elapsed output = Objective_function(guess_init, package, package_moments; nsims = 20000, ret = 1, statcheck = 1)
CSV.write("$dir\\Simulated_data\\simulated_data_base_statcheck.csv", DataFrame(output[2], :auto), header=false) #write CSV file

data_clone = copy(output[2])
kid_pctiles = output[4]
data_clone[:,24] = data_clone[:,24]./data_clone[:,4]
data_clone[:,24] = data_clone[:,24].*data_clone[:,3] #recode

nrows = length(data_clone[:,1]) #useful for later
ncols = length(data_clone[1,:]) #useful for later

ranks = zeros(nrows, 1)
for i = 1:nrows
    ranks[i,1] = findmin(abs.(kid_pctiles .- data_clone[i,24]))[2]
end

data_clone[:,26] = ranks[:,1] #reset child percentile values
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_movecosts_nob_statcheck.csv", DataFrame(data_clone, :auto), header=false) #write CSV file

########################Counterfactdual -- prohibitive moving costs########################
guess_cfact_move = copy(guess_init)
guess_cfact_move[10] = 750 #have to update this one
guess_cfact_move[19] = 750 #have to update this one too
prim_2000_cfact_move, res_2000_cfact_move = Initialize(package, package_moments)
prim_2010_cfact_move, res_2010_cfact_move = Initialize(package, package_moments; p2 = 1)
grids_cfact_move, est_cfact_move = Initialize_grids(prim_2000_cfact_move, guess_cfact_move) #initialize computing grids
V_iterate(prim_2000_cfact_move, est_cfact_move, res_2000_cfact_move, grids_cfact_move) #compute policy and value function.
V_iterate(prim_2010_cfact_move, est_cfact_move, res_2010_cfact_move, grids_cfact_move) #compute policy and value function.
data_simul_cfact_move = Simulate_cfact(prim_2000_cfact_move, prim_2010_cfact_move, est_cfact_move, res_2000_cfact_move, res_2010_cfact_move, grids_cfact_move, package, kid_pctiles; nsims = 20000, statcheck = 1) #simulate CHKS cohort, p25
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_movecosts_statcheck.csv", DataFrame(data_simul_cfact_move, :auto), header=false) #write CSV file

guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112]
prim_2000, res_2000 = Initialize(package, package_moments)
prim_2010, res_2010 = Initialize(package, package_moments; p2 = 1)
grids, est = Initialize_grids(prim_2000, guess_init) #initialize computing grids
V_iterate(prim_2000, est, res_2000, grids) #compute policy and value function.
V_iterate(prim_2010, est, res_2010, grids) #compute policy and value function.

data_simul_decomp_1 = Simulate_decomp_1(prim_2000, prim_2010, est, res_2000, res_2010, res_2000_cfact_move, res_2010_cfact_move, grids, package, kid_pctiles; nsims = 20000, statcheck = 1)
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_movecosts_decomp_1_statcheck.csv", DataFrame(data_simul_decomp_1, :auto), header=false) #write CSV file

data_simul_decomp_2 = Simulate_decomp_2(prim_2000, prim_2010, est, res_2000, res_2010, res_2000_cfact_move, res_2010_cfact_move, grids, package, kid_pctiles; nsims = 20000, statcheck = 1)
CSV.write("$dir\\Simulated_data\\simulated_data_cfact_movecosts_decomp_2_statcheck.csv", DataFrame(data_simul_decomp_2, :auto), header=false) #write CSV file




###################################################################################################################
#######################Amenity Options#######################
###################################################################################################################

################Estimation Procedure

#amenity option 1
function wrapper(guess::Array{Float64}, grad::Vector)
    err = Objective_function(guess, package, package_moments; nsims = 20000, ret = 0, aopt = 1)
    err
end

opt = Opt(:LN_SBPLX, 31)
opt.xtol_rel = 1e-4
opt.ftol_rel = 1e-4
opt.min_objective = wrapper
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112, 
0.008, 0.104]
 #d2shore, cooling
(minf, minx, ret) = NLopt.optimize(opt, guess_init)
println("All finisheed. guess_init = ", round.(minx, digits = 5))
println(minf)


#amenity option 2
function wrapper(guess::Array{Float64}, grad::Vector)
    err = Objective_function(guess, package, package_moments; nsims = 20000, ret = 0, aopt = 2)
    err
end

opt = Opt(:LN_SBPLX, 31)
opt.xtol_rel = 1e-4
opt.ftol_rel = 1e-4
opt.min_objective = wrapper
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112, 
-0.071, -0.175]
 #crime, est
(minf, minx, ret) = NLopt.optimize(opt, guess_init)
println("All finisheed. guess_init = ", round.(minx, digits = 5))
println(minf)


#amenity option 3
function wrapper(guess::Array{Float64}, grad::Vector)
    err = Objective_function(guess, package, package_moments; nsims = 20000, ret = 0, aopt = 3)
    err
end

opt = Opt(:LN_SBPLX, 31)
opt.xtol_rel = 1e-4
opt.ftol_rel = 1e-4
opt.min_objective = wrapper
guess_init = [0.56573, 0.55223, -0.5444, 0.42732, 3.62346, 0.93305, 0.38092, 0.32518, 0.25482,
14.25803, 0.65323, 3.54421, 0.99057, 1.98836,
-1.80509, 0.84117, 0.64976, 1.09643,
15.65022, 1.23573,
-1.57715, -1.99212, 0.11161, -0.25218,  3.55198, 3.58337,
-0.35679, -0.49631, -0.48112, 
-0.159, -0.011]
 #rtw, debt
(minf, minx, ret) = NLopt.optimize(opt, guess_init)
println("All finisheed. guess_init = ", round.(minx, digits = 5))
println(minf)




#############Standard Errors
function standard_errors(guess::Vector{Float64}, perturb::Float64, aopt::Int64)
    num_sim = 20000 #10x normal simulation
    nrep = 1
    num_param = length(guess)
    num_moment = 234
    perturb_mat = zeros(num_moment, num_param)
    se_vec = zeros(num_param)

    #setup
    @everywhere package = Readin_utilities(dir)
    @everywhere package_moments = Readin_moments(dir)
    err_vec = Objective_function(guess, package, package_moments; nsims = num_sim, ret = 1, aopt = aopt)[5] #


    #loop over parameters
    for i = 1:num_param
        guess_perturb = copy(guess)
        guess_perturb[i] += perturb
        err_vec_perturb = Objective_function(guess_perturb, package, package_moments; nsims = num_sim, ret = 1, aopt = aopt)[5] #
        perturb_mat[:, i] = (err_vec.-err_vec_perturb)./perturb
    end

    gradient_mat = perturb_mat
    W = package[14]
    avar_mat_step1 = Matrix{Float64}(gradient_mat' * W * gradient_mat)
    avar_mat_step2 = pinv(avar_mat_step1)
    avar_mat_step3 = (1 + 1/nrep).*avar_mat_step2
    se_vec = sqrt.(diag(avar_mat_step3))
    se_vec
end

#get standard erors
guess_init_1 = [0.56296, 0.55278, -0.55461, 0.42418, 3.62281, 0.93305, 0.38301, 0.32731, 0.25482, 14.24547, 0.66647, 3.53154, 1.00147, 1.98725, -1.80772, 0.84362, 0.64957, 1.04465, 15.68491, 1.23843, -1.57923, -1.96771, 0.16809, -0.23186, 3.5847, 3.58059, -0.35853, -0.50732, -0.49342, 0.01599, 0.06579]
guess_init_2 = [0.56429, 0.552, -0.5444, 0.42648, 3.72291, 0.92609, 0.38361, 0.32494, 0.24917, 14.2395, 0.66105, 3.5443, 0.99106, 1.98632, -1.8085, 0.86847, 0.64529, 1.09427, 15.67903, 1.23865, -1.57989, -1.99776, 0.15644, -0.20326, 3.6538, 3.62802, -0.35929, -0.49018, -0.47918, -0.07001, -0.21876]
guess_init_3 = [0.56034, 0.55966, -0.54823, 0.42885, 3.62693, 0.93471, 0.38346, 0.32406, 0.25477, 14.2513, 0.6605, 3.54631, 0.99766, 1.98712, -1.80644, 0.84352, 0.64932, 1.08866, 15.67067, 1.23638, -1.54773, -1.98025, 0.17566, -0.18034, 3.55193, 3.57348, -0.35166, -0.51547, -0.48708, 0.04458, -0.05267]

se_vec_1 = standard_errors(guess_init_1, 1e-3, 1)
se_vec_2 = standard_errors(guess_init_2, 1e-3, 2)
se_vec_3 = standard_errors(guess_init_3, 1e-3, 3)



println("#############################")
for i = 1:length(guess_init_1)
    println(round(guess_init_1[i], digits=3), " || ", round(se_vec_1[i], digits=3))
end
println("#############################")

println("#############################")
for i = 1:length(guess_init_2)
    println(round(guess_init_2[i], digits=3), " || ", round(se_vec_2[i], digits=3))
end
println("#############################")

println("#############################")
for i = 1:length(guess_init_3)
    println(round(guess_init_3[i], digits=3), " || ", round(se_vec_3[i], digits=3))
end
println("#############################")





##################Main counterfactual exercise

guesses = zeros(3, 31)

guess_init_1 = [0.56296, 0.55278, -0.55461, 0.42418, 3.62281, 0.93305, 0.38301, 0.32731, 0.25482, 14.24547, 0.66647, 3.53154, 1.00147, 1.98725, -1.80772, 0.84362, 0.64957, 1.04465, 15.68491, 1.23843, -1.57923, -1.96771, 0.16809, -0.23186, 3.5847, 3.58059, -0.35853, -0.50732, -0.49342, 0.01599, 0.06579]
guess_init_2 = [0.56429, 0.552, -0.5444, 0.42648, 3.72291, 0.92609, 0.38361, 0.32494, 0.24917, 14.2395, 0.66105, 3.5443, 0.99106, 1.98632, -1.8085, 0.86847, 0.64529, 1.09427, 15.67903, 1.23865, -1.57989, -1.99776, 0.15644, -0.20326, 3.6538, 3.62802, -0.35929, -0.49018, -0.47918, -0.07001, -0.21876]
guess_init_3 = [0.56034, 0.55966, -0.54823, 0.42885, 3.62693, 0.93471, 0.38346, 0.32406, 0.25477, 14.2513, 0.6605, 3.54631, 0.99766, 1.98712, -1.80644, 0.84352, 0.64932, 1.08866, 15.67067, 1.23638, -1.54773, -1.98025, 0.17566, -0.18034, 3.55193, 3.57348, -0.35166, -0.51547, -0.48708, 0.04458, -0.05267]

guesses[1,:] = guess_init_1
guesses[2,:] = guess_init_2
guesses[3,:] = guess_init_3

for i = 3:3
    package = Readin_utilities(dir)
    package_moments = Readin_moments(dir)

    guess_init = Vector(guesses[i, :])
    @elapsed output = Objective_function(guess_init, package, package_moments; nsims = 20000, ret = 1, aopt = i)
    CSV.write("$dir\\Simulated_data\\simulated_data_base_amen_$i.csv", DataFrame(output[2], :auto), header=false) #write CSV file


    data_clone = copy(output[2])
    kid_pctiles = output[4]
    data_clone[:,24] = data_clone[:,24]./data_clone[:,4]
    data_clone[:,24] = data_clone[:,24].*data_clone[:,3] #recode

    nrows = length(data_clone[:,1]) #useful for later
    ncols = length(data_clone[1,:]) #useful for later

    ranks = zeros(nrows, 1)
    for i = 1:nrows
        ranks[i,1] = findmin(abs.(kid_pctiles .- data_clone[i,24]))[2]
    end

    data_clone[:,26] = ranks[:,1] #reset child percentile values
    CSV.write("$dir\\Simulated_data\\simulated_data_cfact_movecosts_nob_amen_$i.csv", DataFrame(data_clone, :auto), header=false) #write CSV file


    guess_cfact_move = copy(guess_init)
    guess_cfact_move[10] = 750 #have to update this one
    guess_cfact_move[19] = 750 #have to update this one too
    prim_2000_cfact_move, res_2000_cfact_move = Initialize(package, package_moments)
    prim_2010_cfact_move, res_2010_cfact_move = Initialize(package, package_moments; p2 = 1)
    grids_cfact_move, est_cfact_move = Initialize_grids(prim_2000_cfact_move, guess_cfact_move) #initialize computing grids
    V_iterate(prim_2000_cfact_move, est_cfact_move, res_2000_cfact_move, grids_cfact_move; aopt = i) #compute policy and value function.
    V_iterate(prim_2010_cfact_move, est_cfact_move, res_2010_cfact_move, grids_cfact_move; aopt = i) #compute policy and value function.
    data_simul_cfact_move = Simulate_cfact(prim_2000_cfact_move, prim_2010_cfact_move, est_cfact_move, res_2000_cfact_move, res_2010_cfact_move, grids_cfact_move, package, kid_pctiles; nsims = 20000, aopt=i) #simulate CHKS cohort, p25
    CSV.write("$dir\\Simulated_data\\simulated_data_cfact_movecosts_amen_$i.csv", DataFrame(data_simul_cfact_move, :auto), header=false) #write CSV file
end


##############################





######
